{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "53ce0ff6",
   "metadata": {},
   "source": [
    "### 2.1. PAO分布\n",
    "\n",
    "PAO分布是指随机变量X的分布，其需要达到的目标是：\n",
    "1. 在偏度为0，峰度为3的情况下，X的分布是正态分布；\n",
    "2. 在偏度为0，峰度为1.8的情况下，X的分布是均匀分布；\n",
    "3. 在标准差一定的情况下，峰度越靠近1.8，分布约接近均匀分布；\n",
    "4. 在标准差一定的情况下，峰度越远离1.8，分布的高度越高；\n",
    "\n",
    "PAO分布的形式化定义为：\n",
    "其概率密度函数为：\n",
    "\n",
    "$$ a = (\\delta - 1.772) * 0.81433 $$\n",
    "$$ z = \\frac {x - \\mu} {\\sigma} $$\n",
    "$$ f(x)=\\frac {1} {\\sqrt{2\\pi}\\sigma} \\exp\\left(-{a{z}^2+b{z}^3+c{z}^4}\\right) $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "09252cf3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.colors as mcolors\n",
    "from scipy import stats, special\n",
    "from scipy.optimize import minimize_scalar\n",
    "from collections.abc import Callable\n",
    "from scipy.optimize import minimize\n",
    "from scipy.optimize import least_squares\n",
    "from scipy.integrate import quad\n",
    "from scipy.stats import johnsonsu\n",
    "from scipy import special as sc\n",
    "# 解决中文乱码问题\n",
    "plt.rcParams['font.sans-serif'] = ['SimHei'] # 设置字体为黑体\n",
    "plt.rcParams['axes.unicode_minus'] = False # 解决负号显示为方块的问题\n",
    "# pao_normal\n",
    "_norm_pdf_C = np.sqrt(2*np.pi)\n",
    "_norm_pdf_logC = np.log(_norm_pdf_C)\n",
    "\n",
    "def _norm_pdf(x):\n",
    "    return np.exp(-x**2/2.0) / _norm_pdf_C\n",
    "\n",
    "def _norm_cdf(x):\n",
    "    return sc.ndtr(x)\n",
    "\n",
    "\n",
    "def _norm_logcdf(x):\n",
    "    return sc.log_ndtr(x)\n",
    "\n",
    "\n",
    "def _norm_ppf(q):\n",
    "    return sc.ndtri(q)\n",
    "\n",
    "\n",
    "def _norm_sf(x):\n",
    "    return _norm_cdf(-x)\n",
    "\n",
    "\n",
    "def _norm_logsf(x):\n",
    "    return _norm_logcdf(-x)\n",
    "\n",
    "\n",
    "def _norm_isf(q):\n",
    "    return -_norm_ppf(q)\n",
    "class paonormal_gen(type(johnsonsu)):\n",
    "    r\"\"\"A Johnson SU continuous random variable.\n",
    "\n",
    "    %(before_notes)s\n",
    "\n",
    "    See Also\n",
    "    --------\n",
    "    johnsonsb\n",
    "\n",
    "    Notes\n",
    "    -----\n",
    "    The probability density function for `johnsonsu` is:\n",
    "\n",
    "    .. math::\n",
    "\n",
    "        f(x, a, b) = \\frac{b}{\\sqrt{x^2 + 1}}\n",
    "                     \\phi(a + b \\log(x + \\sqrt{x^2 + 1}))\n",
    "\n",
    "    where :math:`x`, :math:`a`, and :math:`b` are real scalars; :math:`b > 0`.\n",
    "    :math:`\\phi` is the pdf of the normal distribution.\n",
    "\n",
    "    `johnsonsu` takes :math:`a` and :math:`b` as shape parameters.\n",
    "\n",
    "    The first four central moments are calculated according to the formulas\n",
    "    in [1]_.\n",
    "\n",
    "    %(after_notes)s\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    .. [1] Taylor Enterprises. \"Johnson Family of Distributions\".\n",
    "       https://variation.com/wp-content/distribution_analyzer_help/hs126.htm\n",
    "\n",
    "    %(example)s\n",
    "\n",
    "    \"\"\"\n",
    "    def _pdf(self, x, a, b):\n",
    "        # paonormal.pdf(x, a, b) = b / sqrt(x**2 + 1) *\n",
    "        #                          phi(a + b * log(x + sqrt(x**2 + 1)))\n",
    "        x2 = x*x\n",
    "        trm = _norm_pdf(a + b * np.arcsinh(x))\n",
    "        return b*1.0/np.sqrt(x2+1.0)*trm\n",
    "\n",
    "    def _cdf(self, x, a, b):\n",
    "        return _norm_cdf(a + b * np.arcsinh(x))\n",
    "\n",
    "    def _ppf(self, q, a, b):\n",
    "        return np.sinh((_norm_ppf(q) - a) / b)\n",
    "\n",
    "    def _sf(self, x, a, b):\n",
    "        return _norm_sf(a + b * np.arcsinh(x))\n",
    "\n",
    "    def _isf(self, x, a, b):\n",
    "        return np.sinh((_norm_isf(x) - a) / b)\n",
    "\n",
    "    def _stats(self, a, b, moments='mv'):\n",
    "        # Naive implementation of first and second moment to address gh-18071.\n",
    "        # https://variation.com/wp-content/distribution_analyzer_help/hs126.htm\n",
    "        # Numerical improvements left to future enhancements.\n",
    "        mu, mu2, g1, g2 = None, None, None, None\n",
    "\n",
    "        bn2 = b**-2.\n",
    "        expbn2 = np.exp(bn2)\n",
    "        a_b = a / b\n",
    "\n",
    "        if 'm' in moments:\n",
    "            mu = -expbn2**0.5 * np.sinh(a_b)\n",
    "        if 'v' in moments:\n",
    "            mu2 = 0.5*sc.expm1(bn2)*(expbn2*np.cosh(2*a_b) + 1)\n",
    "        if 's' in moments:\n",
    "            t1 = expbn2**.5 * sc.expm1(bn2)**0.5\n",
    "            t2 = 3*np.sinh(a_b)\n",
    "            t3 = expbn2 * (expbn2 + 2) * np.sinh(3*a_b)\n",
    "            denom = np.sqrt(2) * (1 + expbn2 * np.cosh(2*a_b))**(3/2)\n",
    "            g1 = -t1 * (t2 + t3) / denom\n",
    "        if 'k' in moments:\n",
    "            t1 = 3 + 6*expbn2\n",
    "            t2 = 4*expbn2**2 * (expbn2 + 2) * np.cosh(2*a_b)\n",
    "            t3 = expbn2**2 * np.cosh(4*a_b)\n",
    "            t4 = -3 + 3*expbn2**2 + 2*expbn2**3 + expbn2**4\n",
    "            denom = 2*(1 + expbn2*np.cosh(2*a_b))**2\n",
    "            g2 = (t1 + t2 + t3*t4) / denom - 3\n",
    "        return mu, mu2, g1, g2\n",
    "\n",
    "paonormal = paonormal_gen(name='paonormal')\n",
    "\n",
    "\n",
    "# def pao_normal(mean: float, stddev: float, skew: float, kurt: float = 0) -> None:\n",
    "#     \"\"\"生成jsonson su分布\"\"\"\n",
    "#     return create_dist_by_mvsk(paonormal, mean, stddev, skew, kurt, [0.0, 1.0], ([-np.inf, 1e-6], [0, np.inf]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "b5bbbff0",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\榴弹炮\\AppData\\Local\\Temp\\ipykernel_40436\\4250673353.py:47: UserWarning: No artists with labels found to put in legend.  Note that artists whose label start with an underscore are ignored when legend() is called with no argument.\n",
      "  plt.legend()\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA08AAAIfCAYAAACox+dZAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjEsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvc2/+5QAAAAlwSFlzAAAPYQAAD2EBqD+naQAAcx1JREFUeJzt3Qd0VNXaxvFn0juQEHrovXcQaSqIgqCCImBF7OXq1Wv71Kt4vbZr7woKWAALKjZULKggvffeWyBACiF15lv7HELoBkhypvx/ax0nczITXsxJmGf23u92eTwejwAAAAAAJxV08k8DAAAAAAzCEwAAAAAUAeEJAAAAAIqA8AQAAAAARUB4AgAAAIAiIDwBAAAAQBEQngAAAACgCAhPAAAAAFAEhCcAAAAAKALCEwAg4OzcuVN5eXmH7n/33XeaPn36KX+dCRMmqGfPnlq9enUxVwgA8EYuj8fjcboIAABK04ABA7R79279/vvv1v3LLrtM27dv17Rp007p6wwaNEizZs3SypUrFRoaWkLVAgC8BSNPAABHbdiwQS6X69ARGxurrl27avbs2Uc87osvvrA+/8cffxzzNbZu3aqLL75Y0dHRqlSpkh599FG53e7j/nkrVqzQV199pauvvvrQuccff1wzZszQRx99dNzn7N27V8uWLdPatWutes2xYMECq6bBgwdbf37B+YLDPHbp0qXKyMg44/9HAADvEOJ0AQAAGP/97391/vnnWyNC7777rrp3767ly5erevXq1ucnT5586NaEqwL5+fm68MILlZ6errFjx1qh5f7771dUVJQeeuihY/4cc65+/fq67rrrDp1r2rSpbr31Vt15553q3LmzatasecRzzLS+w8PW4Z566inrOJHffvvN+rsAAHwfI08AAK9Qu3ZttW3bVhdccIHGjx9vjSKNHDny0OdNaAoLC9NPP/10xPO+//57LV68WJ988ok1+nTPPffo9ttv1wsvvHDM6NO4ceOsUacXX3xRISFHvn/49NNPq0KFClYQ27Zt2xGfu+KKK5SZmWmtkzKz3c3oVXh4uD744APrvnnu3LlzrY+PPghOAOA/CE8AAK9jQpIJU2YUyVi/fr318T/+8Q/NmTPHmkZXwKxbio+PV/v27Q+dM00cUlJSrGlzBcxUOjOydP3111sB6WhmuuCkSZOskNSxY0dNnTr10OfMeqbIyEgFBwdbgciMUpkRqoLRqJdeeumIPwsA4J8ITwAAr7Rjxw5Vrlz50KiTCS/33Xefdf/XX3899Lh169apRo0aRzy3YKqfCV1GcnKyevfurcTERA0dOtRav2RGj44+cnJy9N5776lcuXLq1q2b7rrrrmPqevjhh631WGaq3pYtW6zDrMVKS0s7dN8cpn4AgH9hzRMAwKuYNU9vvPGGNm3apMsvv/xQeOrUqZM1ra5ly5bW1D3TMc8wa53MFL/DFdw3nzP+9a9/WeHGTO0zI0YnY0Larl27rKl/h49mGW+99ZY1Rc/o0KHDEZ+74447rKOAWTdVEN4AAP6BkScAgFcwXevMCI4ZHTJrkkaMGGEFFLNuyYw0FawdOueccw41jziRgl04zNczzNop04bcBDAzLc80mTjrrLN0yy23HLE+6YknnlC1atWs8DV69GhdeeWVh77mM888o9tuu0033njjoZCXm5trHVWrVtXHH3986L4Jf2ZNFADAvxCeAABewYST+fPna82aNdaapmHDhlnnTSOGPXv2WO3HTRgyjSAK1kAZcXFxx7QD379//6HPFayhMlP7zPPNyJJhOvk1a9bsiOeZluMmPB3NTBc0Xfqee+45a92VYdY/maYTBY0ngoKCDt0v+BgA4F8ITwAAr2DCjZmSV6dOHSt8FDCjTCYEmWBljpkzZ1rBpaDrnnn8xo0bj/haZsqfYZpOHM+UKVO0b98+axTr6PCUlJR0zOMffPBBa8pfwZorAEBgIjwBALxawXonE6zMYdYhtWrV6tDUPROAzEjV9OnTDz3HBCuzPqpRo0bHfD0zZc80fTDT9o7+/InCU0JCggYOHGh9fKLNdw9XlMcAAHwPcwoAAF7LrE/666+/9Nhjjx1x3myS+/7771tBqFevXlaYMmumXn75ZWs6X0Fjh4I1TwWysrJ01VVXWSNY5usezuzhZJ57vPB0uOzs7BN+zjz/ww8/1BdffGGt3QIA+BdGngAAXsvs4WTahx/dIc+EJzPtbtasWdYUv2+//dYKUKbBg1mXZEaWTIe9w5mRqi5duui7776zpuC1bt3aOm/aij/yyCPWBrumI9/RXfSKEp7MSJNpOGHWSxVszmu+JgDAv7g8BS2JAADwU6Y5xNlnn6369evr7bfftqb/FTD/DDZv3tzqjmfak5t9oE5VxYoV9fzzz1ub5ppue2ZTXQCA/2HaHgDA75m1TXPmzFGtWrWOmcpn7i9evPiMvr4ZjSoYkSI4AYD/YuQJAAAAAIqANU8AAAAAUASEJwAAAAAoAsITAAAAABQB4QkAAAAAiiBguu2ZPTe2bdum2NjYYzotAQAAAAgcHo9H6enpqlKlirVfYFEFTHgywenvdo0HAAAAEDg2b95sbXBeVAETnsyIk7F+/XrFx8c7XQ4cYDau/Omnn3T++eezD0uA4hqAwXUArgFwDWDPnj3W3n8FGaGoAiY8FUzVM/+D4uLinC4HDv2ijIqKsr7//KIMTFwDMLgOwDUArgHk5uZat6e6nIeGEQAAAABQBIQnAAAAACgCwhMAAAAAFEHArHkCAAAA4D+txvPy8pSfn3/czwcHByskJKTYtygiPAEAAADwGTk5Odq+fbsyMzNP+jjTFKRy5coKCwsrtj+b8AQAAADAJ7jdbmvrITOyZDa4NcHo6NElMyplAtauXbusx9arV++UNsI9GcITAAAAAJ+Qk5NjBaikpCRrZOlEIiMjrTb0GzdutJ4TERFRLH8+DSMAAAAA+JSgIowkFddo0xFfs9i/IgAAAAD4IcITAAAAABQB4QkAAAAAioDwBAAAAABFQHgCAAAA4FM8Hk+xPManwtPu3btVq1Ytbdiw4ZSed8UVV+jOO+8ssboAAAAAeJ/Q0FDr9u82yD38MQXP8el9nkxwuuiii045OH3//feaMmWKVq5cWWK1AQAAAPA+wcHBKlu2rJKTk637Zq+n422Sa4KTeYx5rHmOz4enQYMGaciQIZo5c2aRn7N//37ddtttevrpp63/EQAAAAACS6VKlazbggB1IiYvFDzW58PTiBEjrCl7d911V5GfM3z4cGuH4JCQEE2ePFnnnXdeiWx+BQAAAMA7uVwuVa5cWRUqVFBubu5xH2Om6p10xGnHIt8KTyY4nYqNGzfqlVdeUdu2bbVu3Tq9/PLLqlatmr766qvjBqjs7GzrKJCWlmbdmv/BJ/qfDP9W8H3n+x+4uAZgcB2AawBcA/4j+AQBye12W8cRcvbLtexLBc0brdAN80/rz3N5SqINxakU4HJp/fr1qlmz5kkf98QTT2jkyJFatWqVIiIilJ6erho1amj8+PE6//zzj3n8448/bo1UHW3s2LHW3EgAAAAA/i/2wBbVTPlNSXumKTTfbiKxLztI5Z7Zp9TUVMXFxXn/yNOp2rJli3r06GEFJyM2Nlb16tXTmjVrjhueHnroId1zzz1HjDwlJSXpnHPOUUJCQqnWDu9g3l0y0z179uxZrF1X4Du4BmBwHYBrAFwDASAvW64V31ijTEGbZxw67SlXS+5W1ygrqZf0TNNT/rI+E57MFL3ly5cfum+G4Uygqlq16nEfHx4ebh1HMz8g/JAENq4BcA3A4DoA1wC4BvxQylpp7mhpwcdSZop9zhUsNewttb1erlrdFRwUpNCUg5/z9fBkRogiIyOPuZAvv/xya73ThAkT1KFDB7322mvWuwZmNAoAAABAgMrPlVZOkua8L637rfB8XFWpzXVSq6uluMrF8kd5XXhq3ry51QzikksuOeJ8o0aNNG7cOD366KPWuqe6detq4sSJio6OdqxWAAAAAA7Zt1ma94F9ZOw4eNIl1e1hjTKp3vlScPHGHcfD09H9Kk62aW6/fv2sAwAAAEAAcudLa36xR5lW/yh5DnbUi060R5jaXCuVO3kjOp8OTwAAAABwUuk7pfkfSnPHSKmbCs/X7GKPMjW8SAoJU0kjPAEAAADwPh6PtP4Pe5RpxbeSO88+H1FWanmlvZ4psX6plkR4AgAAAOA9MvdIC8ZKc0dJKWsKz1drb48yNblECo10pDTCEwAAAADnR5k2z7JHmZZ+KeVn2+fDYqTmV0hth0qVmjldJeEJAAAAgEOy0qRFn0hzRknJSwvPm6DUdpjU7DIpPFbegvAEAAAAoHRtW2CPMi3+XMrdb58LiZCaXmZPzavaWnK55G0ITwAAAABKXk6mtGSCHZq2zSs8X76BHZhaXCFFlpM3IzwBAAAAKDnJy+1peQvHS9mp9rmgUKlxP3tqXo1OXjnKdDyEJwAAAADFKy9bWva1Pcq06a/C82Vr2M0fWl4lxSTK1xCeAAAAABSPPeukuaOl+R9JmSn2OVew1OBCOzTVPlcKCpKvIjwBAAAAOH35udKqH+xRprW/Fp6PrSK1uVZqdbVUpqr8AeEJAAAAwKlL3SLN+8A+0rcfPOmS6p5nN4Co10sK9q+44V9/GwAAAAAlx51vjy6ZUSYz2uRx2+ejykutr5ZaXyvF15K/IjwBAAAAOLmMZGn+h/Z6pn2bCs/X7GKvZWrYVwoJk78jPAEAAAA4lscjbfjTHmVa/o3kzrPPR5SRWl4ptRkqJdZXICE8AQAAACiUuUdaOM4OTSlrCs9Xa2evZWpyqRQaqUBEeAIAAAACnRll2jLbDkxLvpDys+3zYTFS84H2KFPl5gp0hCcAAAAgUGWlSYs/leaMknYuKTxfsZnU7nqp2eVSeKyTFXoVwhMAAAAQaLYvtEeZFn0m5e63z4VESE0H2FPzqraRXC6nq/Q6hCcAAAAgEORkSku/sEPT1rmF58vXtwNTi0FSZDknK/R6hCcAAADAn+1aaQemBeOk7FT7XFCo1LifHZpqnM0oUxERngAAAAB/k5dttxc3a5k2Ti08X7aGvS9Ty6ukmEQnK/RJhCcAAADAX+xZb29kO/8jKXO3fc4VJDXobYem2udKQUFOV+mzCE8AAACAL8vPk1b9YE/NW/tL4fnYylLra6XW10hlqjpZod8gPAEAAAC+KHWrNO8Dad4YKX174fk659lrmepfIAXzcr848X8TAAAA8BVut7T2V3uUadUkyeO2z0eVl1pdJbW5Voqv7XSVfovwBAAAAHi7jGR7HZNZz7RvY+H5Gp3ttUyN+koh4U5WGBAITwAAAIA38nikDVPtUSbTOc+da5+PKCO1GGKHpsQGTlcZUAhPAAAAgDfJ3CMtHG+HppTVheertrXXMjW5VAqLcrLCgEV4AgAAALxhlGnLHDswLf1Cysuyz4dGS80H2qNMlVs4XWXAIzwBAAAATslOlxZ9am9mu3Nx4fmKTe1RpmaXSxFxTlaIwxCeAAAAgNK2fZE9yrT4Myknwz4XEiE16W+HpmptJZfL6SpxFMITAAAAUBpyMqWlX9qhaeucwvMJ9ezA1GKQFBXvZIX4G4QnAAAAoCTtWmlPy1s4VspKtc8FhdjtxU1oqtmFUSYfQXgCAAAAiltejrTiGzs0bfiz8HzZ6lKb66RWV0sxFZysEKeB8AQAAAAUl70b7I1s530oZe62z7mCpPoX2KNMdc6VgoKdrhKnifAEAAAAnIn8PGnVD/ZaprW/mr7j9vmYSlKba6XW10hlqjldJYoB4QkAAAA4HalbpHkf2Ef69sLztc+R2g2zR5uCQ52sEMWM8AQAAAAUlTtfWvOzvZZp9Y+Sx22fjyovtbrKHmmKr+10lSghhCcAAADg76TvsNcxzRsjpW4uPG865bUdKjW8SAoJd7JClALCEwAAAHA8bre0foq9lmnlJMmdZ5+PKCu1vNLumpdY3+kqUYoITwAAAMDhMnZJCz6W5o6yu+cVSOpojzI1vlgKjXSyQgRieNq9e7fatWun3377TTVr1izy83Jzc9W6dWu99tpr6t69e4nWCAAIMLlZUsoaadcKe2PL9G32C6n9yVJOppSfY7/7bDa4DI+RwmKlqHi7k5Y5ytaQKjaRytWSgoKc/tsAKCqPR9ow1Q5My76W3Ln2+fA4qcUgqc1QqWJjp6tEoIYnE5wuuugibdhwWJovoueee05LliwpkboAAAEm94C0/k9p41Rp41/StvmFU3PORGiUVLGpVOMsqUZnqXpHKSKuOCoGUJwy90gLx9kNIFJWF56v0trel6lpfyks2skK4UUcC0+DBg3SkCFDNHPmzFN63urVq/X888+f0kgVAACHC87PlmvJ59Kq7+2uWbmZRz4gooyU2FBKbCCVqS7FVLCPsBi77XBQqP2udHaGlJMu7d9ttyw2h3nxZUaszNfcMss+pr1ib5JZpZXUoLe9sNx8bZfLqf8FQGDzeFQuY7WCv75NWjZRys+2z4dGS80vt0eZqrR0ukp4IcfC04gRI1SrVi3dddddp/S8m2++WQ8++KAmTZp00sdlZ2dbR4G0tLRDU/7MgcBT8H3n+x+4uAag7QutTlm9lnyqkEVZh0574qrJU6ub3NU7yWNGiExgOpNgY0au9qyXa9s8BW36Sy5z7F0vbZ1rH7/+R574OnI37Ct3syuk8vWK5++HIuF3QQDLSlPQks8UPHeUuu5ecei0p0JTuVtfK3fTy6TwWPsk14dfyz3N76/L4zETPJ3jcrm0fv36Io0kjRo1Sq+//rpmzZql8847T48//vgJ1zyZzw0fPvyY82PHjlVUVFSx1A4A8AEejyqkLVL9nd8oYf+qQ6czwipoa7mO2l62jVIja5b4KFBEzh6rjsqpc5WYvlTBnsKpgXui6mhzQhdtKdtBeSFMDwKKW9nMdaq5+zdV3TtdIe4c61yeK0xby3XQxvLnam9UbUaCA0xmZqY1Cy41NVVxcXH+F5527dqlZs2a6ccff1SLFi2s0HSy8HS8kaekpCRt375dCQkJxf73gG+8wzB58mT17NlToaHs9h2IuAYCjMcj18rvFfzn/+RKttfJeoJClV+/j2bmNVSrS+9UaJhDe7Jkp8u19hfrHXDXmp/l8uTb9YVGyd1soNxtb7Sn9aFE8LsgQORkyLX0CwXPGy3XjkWHTnvKN1Bui6s1OTlB51x4CddAgEpJSVHlypVPOTz5TKvyu+++W8OGDbOCU1GEh4dbx9HMDwg/JIGNawBcAwFg82zpp0ekzTMK1zG0u16ujrfLE1leu7//3gpOjl0HofFSi8vtIyNZWvSp1RbZlbzMeqFnDtXuLnW4RarXi659JYTfBX5qx2K7+YP5uTJrEo3gMKnxJVYDCFf1jnLl5SnP/B7gGghYoaf5ffeZ8GSm28XGxuqNN96w7mdkZFjd+h555BFrDRQAAErfKf34f5JpBmGEREpn3W4fpp24N65jMI0oOt1h12jaJM98W1r5vbRuin2UbyB1uVdqOkAK9pl/toHSZbYRWPql3WZ8y+zC8/F17H2ZWgyRopl5hDPndb+FzfS6yMjIY9Kgmdp3dLc+Mxp1wQUXlHKFAACv43ZbjSD082NSVqqZFC61vFI692Eprop8gllvUauLfezdKM0eIc39QNq9UvryJmnKU1Lne6QWg6WQMKerBbxD8go7MJlW49bPvuw92Br1tTvm1erKWib4d3hq3ry5Xn75ZV1yySVHnD96TVRERIQqVaqksmXLlnKFAACvsneD9OWt0qa/7PuVW0p9X/HtNsPlakjnPyl1vU+aNUKa/ob99/zmH9Lvz0nd7rfDISNRCNSNrJd/I815v/Dn3ihb3Q5Mra6yR3SBEuD4b92j+1UUddPcKVOmlFBFAACfYP79WDhe+v4+e12DWdd07iNS+5v8J1SY/aa6/kvqeKs0d7Q07VUpbYsdov56VTrnYXsdB2uiECijTGaE2YwyHdhrn3MFSw0utKfm1T6XnwWUOD/51wUAEFCy0uwAYdY4GEkdpf7vSOX8dAP1sGh7TVTbYdKc96Q/X5BS1kifD5UqvSSd92+pbg+mJ8E/1zKZTWzNmwcFDWCMuGpS62uk1lf7ztRc+AXCEwDAt+xaKY2/UkpZba9t6P6gdPY//We06WRCI+wQZV40Tn9T+us1ybRg/vgyqWYXe6qfL09XBArsWHJwlOkTKTv1yFGmNtdJdcwoU7DTVSIABcC/NAAAv7H0K2ni7db+LYqtIg38QEpqp4ATHit1f0Bqd4M07SVp5rvShj+ld7vbDSXOe5R34+F7cvZLS76wR5m2zjlyLVPra+21TLGVnKwQIDwBAHykm94vw6VpL9v3zSjLZaOkmEQFNNN62Yw2mXVevzwhLf5MWjhWWvaV1Okf0tn/sKf8Ad5s2wJ7lGnRZ4X7MplR5YZ97FGmWt1ZywSvQXgCAHi33APSFzdJy7+275tQcN5jgTFNr6jMO/MDRtqb6pp9rjbPlH5/xn5Beu6j9mgULz7hbesWzX5sc8dI2xcUno+vbY8ytRxCxzx4Jf7lAQB4r/27pXGDpS2zpOAw6eI3peaXO12V96rWVrr+R3uB/eR/S/s2ShNvszfe7fWUvYcU4GSHzK3zpHmjpcUTpNz99nnzs23ty3SdVKMzQR9ejfAEAPBOKWvtRgh71kkRZaVBY6WaZztdlfczHfeaXGIvrJ/5jvTH83ZTiTEXSQ0vkno+ISXUcbpKBBKzee2iT+1Rpp2LC88n1LMDkxkZNVNQAR9AeAIAeGenrQ8uljJ321PSrpwgJdZ3uirfEhJur3kym+lOedreUHTFt9KqH6UON9sb8Eay0TxKcJ3ixqnS/I/skdC8LPt8cLgd7k1oqn4W7fXhcwhPAADvsmWu9FF/KWufVKm5dNUE1j6cCfOOfp/n7c58Pz0srflZmv66vdHoOf8ntb6O9WMoPqlbpAXjpAUfSXs3FJ5PbGQHpuYDpah4JysEzgi/LQEA3mPDNGnsQLsVeVIH6crPpIgyTlflHyo0tIPo6snSjw9Lu1dK390rzRop9XrS3mQXOB152dLK7+1RpjW/mMVN9vmwWKnZAKnV1VLVNowywS8QngAA3sG86DKb3+YdkGp1lQaNk8JjnK7K/9TrKdXubu+l89tT0q7l0kcDpHrn223PExs4XSF8aXqtCUyLPpEO7Ck8b5o+tL5aatRPCotyskKg2BGeAADOWzdFGj/EXhdRr5c0cIwUGul0Vf4rOFRqf6PU7DLp9/9Js96RVv9kB9h2w6TuDzG1Csd3YJ/dYnzeh0e2GDebVpv24uagIQn8GOEJAOCsjX/Z7chNcKp/oTTwAykkzOmqAkNkOemCp6S210uTH7WnXs161x5J6PagvU6K7wXc+dL6P6QFH0vLvyls/hAUKjXsbU/Lq3OuFBTsdKVAiSM8AQCcs3m29PHlUm6mvebGjDjxYr30la8rDR5njwCa9VA7l0g/PiTNNuuh/ivVv4D1KoEoebndWMS0GU/fXni+QmM7MDW/ghbjCDiEJwCAM7YtsNfamOYQZo3TFR/Z7bXhHLMW6uY/pPkfSr8+Ke1ZK40bJNXqZq+Hqtzc6QpR0jKSpcWfS4vGS9sXFp43e601Nc0frpKqtCJMI2ARngAAzmyAa4JTdqpUvZM0eDxrnLyFmXplWko36S/9+YI0401p/e/SO13sc+c8bI9UwX/kHrCnbC78xG5l78kvnJZXv5fUYpDdUIQ3NwDCEwCglKXvlD681N4A1+zjNOQTKSza6apwtIg4qedwqe1Q6ZcnpCUTpKVf2BuemqYA3R6QyiY5XSXOZBPbTdPtEaalX0nZaYWfq9rWDkxmpInGIcARCE8AgNKTlSZ9PEDat1EqV9Ped8i8SIf3Mt+ny96XOv/Tnsq36gd7Wp9pKtF2mNTlHjYx9hUej7Rt/sEg/KWUtrXwc2WqSy2usNcxla/nZJWAVyM8AQBKbyPNT66UdiyWohOlq77gRbcvqdTMHiXcNNMeido4VZr5lr1flJnm1+lOqUxVp6vE8excZgcmc+xdX3g+PE5q3E9qMdiePhsU5GSVgE8gPAEASmeK0Jc32+2Ow2KkKz9nLxhfVb2DdN230rrf7JGorXPtEGU687W6Ujr7bim+ltNVwqwrNNMsl3whJS8rPB8SKTW40J6SZzpchkY4WSXgcwhPAICS98vj9jQhswDddNWr0tLpinAmTKc1s69P7XOktb/ajSU2TrNHoczmqWbz3bPuoDtfaU/J27XC3ofJHDsWFX7O/NzV62kHJtN2PjzGyUoBn0Z4AgCUrHkfSNNesT++5E2pzjlOV4TiDFF1z7MPs9nxH89La3+x10OZo0ZnqeOt9kgHG6iWzIjutnmFgcm0li/gCpZqd7MDU8M+9obIAM4Y4QkAUHLW/yl9+0/7Y9OdrflApytCSanRSbr6C2nrPGn663ZXPrMuyhxla0jtb7L3CIos63Slvi0vR9r0l7T8W2nFd1L6tsLPBYfZo4GN+tqBNbq8k5UCfonwBAAoGbvXSJ9cJbnz7He/uz/kdEUoDVVb2935Urfa66DmjrK7K/70sPTrf6RG/aTWV9ujUjQoKJq07dKaydKqH6V1U+yNpQuYNYRmD6ZGF0l1e9K9EihhhCcAQPHL3CONHShl7ZOqtZMuftOe4oXAYTrv9XhM6nqftPhTaea7UvJS+2NzmBboZiSq6WU0mDje6JJpxGEC0+qf7A6Vh4uuINU3gamfVKsbTR+AUkR4AgAUr/xc6dNr7PUXZu+YQWN5cRfIwqLsVuatr7XX55iGElbL7A12tz5zVGklNekvNblEKltdAcedL21fYE9zNR0pzea1uZmHPcAlVW1jjzCZxg+VWzJqBziE8AQAKF4/PSJt+NOeTmT2BWIvJxiugwHAHL2estdELRpvhwWzcas5Jj8qVWldGBJMqPLHRhNmzzMzmrR5prRhqrRhmpSdeuRjohKk2t2ler3shhysXwK8AuEJAFB8FoyTZr5tf9z/XaliY6crgreORrUcbB8Zu6TlX9ut7E2QMKNT5vj9GTtA1DnP7hpX/SwpvrbvTf80LcTTtkpbZkubZ9u32xdK+dlHPs5sWFuzs1Srq1Szi1ShMaNLgBciPAEAiofpsvbNXYWd9Ux7ZODvxCRK7YbZR/rOwnU+a6dImSmFa6SM6ESpekcpqaO9h1TFplJUvLxGbpa0e5W0c4m0wxyL7I8P7D32sZHxUlJ7++9jAlOlFlIwL8sAb8dPKQDgzJnRA9NZz7ybXv9CqduDTlcEXxRb0W4iYQ6zdm7zLGnNz/YaINNAYf+uwj2NCsRVtUNUhYZ2EwrTFt3clkmSQsKKt778PLuG9O32sWedlLLWvjVH6hYz1HTs88yeSxWb2GHJNFAxhy+OogEgPAEAzpB5kfvZtfbUpIR6Uv93mG6EMxccKtU82z4KRnVMUwUTpLbMsdcMmRbo5rozx+ofj/oCLntUKqq8vV7ITAE090MiFRQcpgbbNynorzVSSKjdTt+TbzduMB/n7Jey0ux1SNnp9siRGRXbnyx53CevO6KsHeYqNS28TWxE0xTATxCeAABn5seHpY3TpLBYu7NeRBmnK4I/MuHDTHEzR4GsVGnnMntq3O7VdpgyXfz2bpTyDtjT/syxe+URX8q0oGhoPtjx5anXYUaRYitJMRWlcjWk+DpSQp3CWxPSGFEC/BbhCQBw+hZ+Is16p7BBRGJ9pytCIDFBvcZZ9nF0kwYzvS4jWcrcLe3fbe89ZkaQ8g4oPydTm9atVo0qFWSNkQaF2KOl1m2IFBZtN3AIj7X/DHOYwBRb2Q5H/tgBEECREJ4AAKcnebn07d2HNYjo7XRFgM2M/JgW+Sdok+/OzdWi779Xtd69FRQaWurlAfBdTEoHAJy67Azp02vtjTzNXjQmPAEA4OcITwCAU2OmRJkRJ7OOxExj6j+SaUwAgIBAeAIAnJo570uLP7MXzl82yt6nBwCAAEB4AgAU3bYF0g8H93Dq8dixC/UBAPBjhCcAQNEc2Cd9eo2UnyM16C11+ofTFQEAUKoITwCAoq1zmni7vY9O2erSJW+ylw0AIOAQngAAf2/669KKb6XgMGngB1JkOacrAgAgsMLT7t27VatWLW3YsKFIj3/33XdVuXJlhYaGqlu3btq+fXuJ1wgAAW/zLGnyY/bHFzwtVWnldEUAAARWeDLB6aKLLipycJo6daoeffRRffjhh1q/fr08Ho/+9a9/lXidAKBAX+f0+TDJky81HSC1HeZ0RQAABF54GjRokIYMGVLkx69evVrvvPOOevTooWrVqmno0KGaP39+idYIAAHN2s/pn1LqJqlsDemil1nnBAAIaCFO/cEjRoywpuzdddddRXq8CUuHW7lyperVq1dC1QEANP8jaekXUlCIdNn7UkSc0xUBABCY4ckEp9O1Z88eaxRq7NixJ3xMdna2dRRIS0uzbnNzc60Dgafg+873P3BxDZyC3asVMul+mXGm/G7/J3fFFuZ/nPwB1wG4BsA1gNzT/N67PGbxkINcLpe1hqlmzZpFfs7gwYOtMPTdd9+d8DGPP/64hg8ffsx5E7iioqJOu14A8HdB7hx1XfWEyhzYpF0xjfVX3fslF81ZAQD+IzMz01pClJqaqri4OP8NT2PGjNEDDzyghQsXqmLFiqc08pSUlGR16EtISCiW2uF77zBMnjxZPXv2tDo2IvBwDRRN0E8PK3j2O/JEJSjvhilSbGX5E64DcA2AawApKSlWF+9TDU+OTds7HXPmzNGdd96pr7/++qTByQgPD7eOo5kfEH5IAhvXALgGTmLVj9Lsd6wPXZe8pdD46vJXXAfgGgDXQOAKPc3vu9fNwzAjRMebg5icnKy+ffvq/vvvV9u2bZWRkWEdAIBikr5D+upW++MOt0r1ezldEQAAXsXrwlPz5s2Pu5Zp3Lhx2rFjh7XXU2xs7KEDAFAM3G7pi5ukzBSpUjOp57FrRgEACHSOT9s7esnViTbNNS3Ni9rWHABwiv56RVr/uxQaJQ14Xwo5dtozAACBzutGngAApWzLHOnXJ+2PL3xWSqzvdEUAAHglwhMABLKsVOnz6yV3ntTkUqnV1U5XBACA1yI8AUCgMtOmv71H2rdRKlNduuhls3+E01UBAOC1CE8AEKgWjpOWfC65gqUBI6XIsk5XBACAVyM8AUAg2r1G+u5f9sfnPCRV7+B0RQAAeD3CEwAEmrwcacL1Uu5+qWYXqfM9TlcEAIBPIDwBQKD5Zbi0faEUWU7q/64UFOx0RQAA+ATCEwAEktU/S9Nftz+++A0prorTFQEA4DMITwAQKNJ3Sl/dYn/c7kapYR+nKwIAwKcQngAgELjddnDav0uq0EQ6/z9OVwQAgM8hPAFAIDBT9db+KoVESpe9L4VGOl0RAAA+h/AEAP5u6zy7SYRxwdNShYZOVwQAgE8iPAGAP8tOlyYMk9x5UqO+UpvrnK4IAACfRXgCAH/2/X3SnnVSXDWp76uSy+V0RQAA+CzCEwD4q4WfSAvHSa4gacAIKSre6YoAAPBphCcA8Ecpa6Xv7rE/7vaAVKOT0xUBAODzCE8A4G/ycux1TjkZUo2zpa73OV0RAAB+gfAEAP7GdNbbNl+KKCv1f1cKCna6IgAA/ALhCQD8yeqf7T2djEvelMpUc7oiAAD8BuEJAPxF+k7pq1vsj9vdKDXs43RFAAD4FcITAPgDt1v68mZp/y6pQhPp/CedrggAAL9DeAIAf/DXq9K636SQSOnyUVJohNMVAQDgdwhPAODrtsyVfv2P/fGFz0qJDZyuCAAAv0R4AgBflpUmTbhecudJTS6VWl/jdEUAAPgtwhMA+CqPR/r2n9LeDVKZ6tJFL0sul9NVAQDgtwhPAOCrFoyVlnwuuYKly96TIss6XREAAH6N8AQAvmj3aun7++yPz31YSmrvdEUAAPg9whMA+JrcA9Jn10m5+6VaXaWz/+l0RQAABATCEwD4mkkPSDuXSNEVpP4jpCB+lQMAUBr4FxcAfMmiT6V5YyS5pAEjpNhKTlcEAEDAIDwBgK/YtUr65m77424PSLW7O10RAAABhfAEAL4gJ1P67NqD65y6Sd3ud7oiAAACDuEJAHzBpPuk5GVSTEVpwEgpKNjpigAACDiEJwDwdgvGSfM/klxBdnCKqeB0RQAABCTCEwB4s+QV0nf32B93f8huTQ4AABxBeAIAb5Wz/+A6p0yp9jlSl3udrggAgIBGeAIAb+TxSN/+U9q1QoqpdHA/J9Y5AQDgJMITAHij2SOlRZ9IrmDpsvelmESnKwIAIOARngDA22yaKf3woP3x+f+Rap7tdEUAAIDwBABeJn2nvc7JnSc1uVTqeJvTFQEAgIMITwDgLfJzpc+HSunbpfINpH6vSy6X01UBAICDCE8A4C1+flzaOE0Ki5Wu+EgKj3G6IgAA4C3haffu3apVq5Y2bNhQpMf//vvvatSokcqXL68XX3yxxOsDgFKz9Etp+uv2x5e8KSXWd7oiAADgLeHJBKeLLrqoyMFp165d6tevnwYPHqzp06fr448/1m+//VbidQJAidu1Uvrqdvvjs++SGvdzuiIAAOBN4WnQoEEaMmRIkR9vwlKVKlX06KOPql69evr3v/+t9957r0RrBIASd2CvNG6wlLtfqtlFOvffTlcEAABOIEQOGTFihDVl76677irS4xcuXKhzzjlHroOLp9u3b68HHzzYyvc4srOzraNAWlqadZubm2sdCDwF33e+/4HL664Bd56CP7teQXvWyhNXTXmXvCu5PZLbS+rzU153HaDUcQ2AawC5p/m9dyw8meB0Kkz4ady48aH7cXFx2rZt2wkf//TTT2v48OHHnDdT/aKiok6xWviTyZMnO10CHOYt10CTLWNVd9evygsK05+Vb1ba77OdLimgeMt1AOdwDYBrIHBlZmb6Vng6VSEhIQoPDz90PyIi4qR/6Yceekj33HPPEeErKSnJGr1KSEgo8Xrhne8wmF+SPXv2VGhoqNPlIMCvAdei8QqZ/4N955K31LnRxY7WE0i86TqAM7gGwDWAlJQU/w5P8fHxVtOIAunp6QoLCzvh403QOjxsFTA/IPyQBDauATh+DWyeLX1/8M2drvcrpPllztUSwBy/DuA4rgFwDQSu0NP8vvvMPk/t2rWzuuwVmD9/vqpWrepoTQBwytK2SZ9cKeXnSA0vkro/5HRFAADAV8OTmV53vAVcpk35tGnT9PPPP1uff+6559SrVy9HagSA05J7QBo/RMrYKVVoLF36jhTkdb+GAQDACXjdv9rNmzfXd999d8x5szHuSy+9pN69e6tixYpauXKlHnnkEUdqBIBT5vFIX98pbZsvRcZLg8dJ4TFOVwUAAHxpzZPHvKA4zMk2zb3lllus0aYVK1aoS5cuionhhQcAHzHlGWnxZ5IrWBo4RipX0+mKAACAr4Wn02lxfqptzgHAUQvGSb8/Y3980UtSra5OVwQAAPxh2h4A+JX1f9rT9YzO/5TaXOt0RQAA4DQRngCgpOxaZXfWc+dKjS+Rzv230xUBAIAzQHgCgJKwf7c09nIpK1Wq1l669G066wEA4OP4lxwASqIl+bjB0t4NdmMI01kvNNLpqgAAwBkiPAFAccrPkybcIG2ZJUWUkYZ8JkWXd7oqAABQDAhPAFBczNYL390jrfhWCg6XBo2VEus7XRUAACgmhCcAKC6/PSXNGyO5gqQBI6WanZ2uCAAAFCPCEwAUh1kjpD+esz/u84LUuJ/TFQEAgGJGeAKAM7X0K+n7++yPuz8ktb3e6YoAAEAJIDwBwJlY/4f0xY1mwZMdmro94HRFAACghBCeAOB0bZ4tjR0k5edIjfpKvZ+XXC6nqwIAACWE8AQAp2P7QumjAVLufql2d6n/SCko2OmqAABACSI8AcCpSl4ufXiplJ0qVT/LbkkeGuF0VQAAoIQRngDgVKSslT64WMpMkaq0loZ8KoVFO10VAAAoBYQnACiqfZukMf2kjJ1SxabSVROkiDinqwIAAKWE8AQARZG61Q5OaVuk8vWlq7+SouKdrgoAAJQiwhMAFGXEaXRvae96qVxN6ZqJUkyi01UBAIBSFlLafyAA+JS9G6TRfaXUTVK5WtJ130pxVZyuCgAAOIDwBAAnsmedHZzMVL34OgQnAAACHOEJAI5n9xppTF8pfZu9xunab6TYSk5XBQAAHER4AoCj7Vp5sKveDimxoR2cYio4XRUAAHAY4QkADrd1nvTRAOnAHqlCE5pDAACAQwhPAFBg3e/S+CFSToZUpZV05QQpOsHpqgAAgJcgPAGAsexracIwKT9HqtVVGjRWCo91uioAAOBF2OcJAOZ9IH12rR2cGvWVhnxGcAIAAMcgPAEIXB6PNPVl6es7JY9banW1dPkYKTTC6coAAIAXYtoegMCUnyf98IA0e6R9/+y7pB7DJZfL6coAAICXIjwBCDymIcRXN0urf5Tkks7/j9TpTqerAgAAXo7wBCCgROTuVcgHfaWdi6WQCKn/u1Lji50uCwAA+ADCE4DAkbxMXVcOlyt3jxRVXho8Xkpq53RVAADARxCeAASG1ZMV8tl1Cs3NkCehrlxXfi7F13K6KgAA4EPotgcgADrqvSR9fLlcORnaHdNQeddOIjgBAIBTxsgTAP+Vkyl9fYe0ZIJ1N7/VNfrL3V0XRpZzujIAAOCDGHkC4J/2bZLe72UHp6AQqc+Lcvd+UR7zMQAAwGngVQQA/7NhqvTptVLmbrsxxMAPpJpnS7m5TlcGAAB8GOEJgP9wu6W/XpV+eULy5EuVW0hXfCyVTXK6MgAAEIjhadOmTZowYYIWLlyo5ORkxcTEqHr16urdu7e6d++uoCBmAgJwQOYe6atbpVU/2PebXyFd9LIUFuV0ZQAAINDC065du/Tggw9q0aJFGjRokG644QZVqFBBGRkZWrt2rUaPHq17771X//vf/9SjR4+SrRoADrdlrvTZtVLqZik4XOr9nNT6WsnlcroyAAAQaOFp1qxZuu666/R///d/eu+99475fOvWrXX55Zdr9erVVqiaMWOGHnnkkZKoFwCObEM+613px4cld64UX1u6fIxUubnTlQEAgEANT2Z0yUzVa9So0UkfV69ePf3888/66KOPiqs+ADi+/buliXdIqybZ9xv1ky5+XYoo43RlAAAgkMPTueeeW+QvmJeXp6FDh55JTQBwcmt+lr68VdqfLAWHST3/I3W4mWl6AACgRJ12d4cuXbpo8eLFR5z77rvv1KRJkyI9f8mSJWrXrp3KlSun++67Tx4z/eYkzOdvvfVWxcfHq2zZstY0wgMHDpxu+QB8UW6W9MND0kcD7OCU2FC68Tep4y0EJwAA4L3hyYxGde3aVXfffbeWL1+uyy67TDfffLOefPLJv31udna2+vbtqzZt2mjOnDlatmyZ1XDiZD788EOtXLlS8+fP159//qmlS5fq6aefPt3yAfiancukEedKM96077e/SbppilSpqdOVAQCAAHHa4Wn48OFatWqVpk+frqZNm2r//v3W/SFDhvztcydNmqTU1FS9+OKLqlOnjp566qnjNqI4ummFCWg1atRQs2bNdMkll2jNmjWnWz4AX5GfK/3xP+mdrlLyUik6URrymdT7f1JopNPVAQCAAHLam+SaMPPAAw8oMzPTCkHPP/+8hg0bZrUqr1at2kmfa/aI6tixo6Ki7P1Xmjdvbo0+nYyZDmhGnwYMGKCsrCyNHz9e99xzz0lHt8xRIC0tzbrNzc21DgSegu87338fkrxMId/cIdeORdZdd71eyu/9khRTwXwjT/nLcQ3A4DoA1wC4BpB7mt97l+fvFhudgFl39Pjjj+vOO+9UcHCw0tPTrVbmY8aMORRUTsTsB2UC0BtvvHHoXGJiojVyZdZAnegvaKb5FayzMtP+vvrqqxNuymtqM6NjRxs7duyh0AbAO7k8eaq38zs12PGVgjz5ygmO1uJqV2tLubNY2wQAAM6YGQAyM+bMbLi4uLiSD087duxQpUqVjjk/e/ZsqxHEyZgRKxOGzIhVgaSkJGt/qKpVqx73OWZk6+uvv7baoLtcLmt9lWmd/sILLxR55Mn8Gdu3b1dCQsIp/E3hL8w1N3nyZPXs2VOhoaFOl4MT2blUId/eedho0wXKv/B5KfbY3zenimsABtcBuAbANYCUlBRVrlz5lMNTSFE745mRoYoVKx46d7zgZJjg9Msvv+i888474dczHfPM1zycGbkKCws74XM+/vhjPfHEE6pevbp13zSL6Nat2wnDU3h4uHUczfyA8EMS2LgGvFTOfun3Z6Xpb0juPCmirLWuKajZ5Qoq5tEmrgEYXAfgGgDXQOAKPc3ve5EaRqxfv17nnHOOpk6detLH7du3T4MGDbKmxp2MCVim0cThX9+MEplQdSJut1vJyclHjHzl5+cXpXwA3m7VT9KbHaVpr9jBqeFF0u2zpOYDmaYHAAC8RpFGnsz6ItMV77bbbrMCzlVXXaUOHTpYI1EZGRlau3atJk6cqE8//dRq4nDTTTed9OuZFudmGt2oUaOsDXVNt70ePXpYa6dMAIuNjbU+PnpfqWeeecY6n5OTo2effVb9+vU7s789AGelbZd+eEBaNtG+H1fN7qLXsLfTlQEAAJx+t73GjRtrypQp+u233zRhwgS9/PLL1khQTEyM1T78/PPPt0aTTtTw4Yg/NCREI0eO1ODBg60Nck3TB/O1DfN8s5dTy5Ytj3iO2T/KBK7777/fmuLXq1cvvfLKK0UtH4A3yc+T5rwn/fIfKSddcgVLHW+Vuj8khcc4XR0AAEDxtCo30/fMcabMqJEZsZo7d67VtrygicOJ+leY7n4ffPDBGf+5ABy2/g9p0oP2nk1G1bZS35elSs2crgwAAKBk9nkqDqbpRJ8+fZwsAUBp2btB+ukRafk39v3IctK5j0hthkpBR07TBQAA8PnwZEaJTLc703kPAIrcRW/qS9K0V6X8bHuKXrth9hS9qBM3iQEAAPA2Req2V8CMEu3Zs+fQGqj9+/eXVF0AfJ3bLS36THq9nfTH/+zgVKurdMtUuykEwQkAAPjzyJNp1FDQEGLFihVW+3AAOMa6KdLkx6TtC+z7ZatL5/9XatSX1uMAACAwwpNpFGFGn5o1sxd233rrrcfd2Pb9998vvgoB+I4di+3QtPYX+35YrNT5LumsO6XQCKerAwAAKL3wNH78eI0ePVopKSlyuVxKSkpSRAQviICAt2+T9Ot/pUWfmJ6ZUlCo1PZ6qdv9UnR5p6sDAAAo/fBk9nS64447rI+HDx+uhx56SHFxccVTCQDfk7lH+vMFada7Un6Ofa5Jf+m8R6X42k5XBwAA4B2tyseNG6fo6OjirQaAbziwT5rxpjTjLSk7zT5nmkH0GC5Vbe10dQAAAN4Vnq644orirQSA98tKk2a+I01/TcpKtc9VbCb1eFyqex7NIAAAgF9zdJNcAD4iO8OemvfXq9KBvfa5xIb2Xk2N+klBp7TrAQAAgE8iPAE4sZxMafZIadrLUmaKfS6hntT9QanJpVJQsNMVAgAAlBrCE4BjZadLc0ZJf70m7U+2z5WrZYemZpcTmgAAQEAiPAE4snuemZ5nGkFk7Svc4LbbA1LzQVIwvzIAAEDg4pUQACl9pzT9dWnO+1JOhn0uvo7U+Z9S8yukkGM3wwYAAAg0hCcg0De3nfaKNO9DKT/bPlexqdTlHqnxJUzPAwAAOAzhCQhEu1ZJU1+SFn8qufPsc9XaSV3+JdXvRctxAACA4yA8AYHC45E2/mU3gVg1qfB8rW5S139JNbsQmgAAAE6C8AT4u/w8acU30rRXpW3zDp50SQ1629PzqrV1uEAAAADfQHgC/FXOfmn+x3YjiH0b7XPB4VLLIdJZt0vl6zldIQAAgE8hPAH+2DnPtBs3m9sWtBuPjJfa3yi1u1GKSXS6QgAAAJ9EeAL8RfIKacYb0sLxUn5O4ca2ZpSp5ZVSWJTTFQIAAPg0whPgy9xuac1ke1Pbdb8Vnq/aVjr7H1LDi2g3DgAAUEwIT4Avyk6XFoyTZr4t7Vlrn3MF2U0gzrpDqt6RznkAAADFjPAE+JK9G6RZI6R5H0jZafa58DJS66vtNU3lajpdIQAAgN8iPAE+sT/TNHtq3srvJY/bPp9QV+pwi9RisBQe43SVAAAAfo/wBHir3CxpyefSjLelnYsLz9c5V+pwq1S3hxQU5GSFAAAAAYXwBHibtO3SnPftI3O3fS40SmoxSGp/s1ShodMVAgAABCTCE+A1U/P+svdnWvGt5M6zz8dVs9cytb5Giop3ukoAAICARngCnJSdIS36xN7QNnlZ4fnqnaQON9utxoP5MQUAAPAGvCoDnLBrlR2YFo4r7JpnpuY1v0Jqd4NUqanTFQIAAOAohCegtOTnSat+kGaPkNZNKTxvuuaZwGS65kWWdbJCAAAAnAThCShpGbukeWOkOaOktC2FG9rWv8Bez1SrO13zAAAAfADhCSipBhBb5tijTEu/lPJz7PNRCXbzh7bXS2WrO10lAAAATgHhCShOuQekJROkWSOk7QsKz1dtI7W7UWpyqRQa4WSFAAAAOE2EJ6A47FkvzXlPmv+RdGCvfS44XGo6QGp/gx2eAAAA4NMIT8Dpcrultb/Yo0yrfzJz9ezzZjpe22FSq6ul6ASnqwQAAEAxITwBp9MAYv6H0txR0r5NhefrnGc3gKh3vhQU7GSFAAAAKAGEJ6CoDSA2TpPmvC8t+1py59rnI8pILa+0W40n1HG6SgAAAJQgwhNwMmb90sLxdmjavarwfNW2dsc80wAiLMrJCgEAAFBKCE/A8UaZts6zA5PpnJd3wD4fGi01v9wOTZVbOF0lAAAAShnhCSiQnSEt/swOTTsWFZ6v0ERqd73UbKAUEedkhQAAAHBQkFN/8JIlS9SuXTuVK1dO9913nzzm3f4icLvd6tSpk1544YUSrxEBYudS6bt7pRcaSt/ebQcn02a8+SBp2GTp1mn2miaCEwAAQEBzJDxlZ2erb9++atOmjebMmaNly5Zp9OjRRXru22+/rdTUVP3jH/8o8Trhx3KzpIWfSO/1kt7qJM0eKeWkS/F1pPP/K927Qur/jpTUXnK5nK4WAAAAgTptb9KkSVYAevHFFxUVFaWnnnpKt99+u4YOHXrS523btk3/93//py+++EKhoaGlVi/8R3TWDgX9/G9p0bjCzWyDQqSGfey1TDW7SkGODcgCAADAizkSnhYuXKiOHTtawclo3ry5Nfr0d+6++27VqFFDmzdv1l9//WVN3zvZ6JY5CqSlpVm3ubm51oEAkp8r1+ofFDRnlHps/OPQaU9cNblbXSN3iyFSbKWDj823D/ilgp99fgcENq4DcA2AawC5p/m9dyQ8mSBTq1atQ/ddLpeCg4O1d+9eaw3U8UyfPl2fffaZevfurbVr1+rJJ59Ur1699Prrrx/38U8//bSGDx9+zPnffvvtUGiDf4vM2a0aKb+rxu7fFZG3zzrnkUs745prQ/lztTOuhZQWJP05z+lSUcomT57sdAnwAlwH4BoA10DgyszMPK3nuTxF7dRQjB544AEr7ZlpewWSkpI0Y8YMVa1a9bjPuf76663RKROiTNgyo09mFGr58uVq0KBBkUaezJ+xfft2JSQklNDfDI5z58u17lcFzRst15rJcnnc1mlPdAXlNRuk39JqqPNFQ5j2GaDM7x3zD2XPnj25BgIY1wG4BsA1gJSUFFWuXNlaShQXF+fdI0/x8fFWt73DpaenKyws7ITP2bJlizXqZIKTYYJQYmKiNQp1vPAUHh5uHUczPyD8kPihjGRp/ofS3NHSvk2F52t2kdoNk6tBH8nj0oHvv+caANcALFwH4BoA10DgCj3N77sj4cm0KB8xYsSh++vXr7dGiUyoOpFq1arpwIEDha+VMzK0Z8+eE45UIQCYQdMNU6U570nLv5XcB+euRpSVWl4ptblOSqxf+HjmNQMAAOAMOBKeunbtak2jGzVqlNVhz3Tb69Gjh7Xuad++fYqNjbU+PtzgwYOtwzyubt26evTRR9WwYUOr2QQCjOmSt2CcvZltyurC89XaSW2HSU0ukUIjnawQAAAAfsiR8BQSEqKRI0daYchskBsUFKQpU6ZYnzMNI+bPn6+WLVse8RwzJ/XZZ5/Vrbfeaq13Mp///PPPD03jQwCMMm2dawemJROkvCz7fFiM1Hyg1GaoVJkgDQAAAD8LT0a/fv2s9Upz58612pYXNHE4Wf+KYcOGWQcCSHa6tPgzOzTtWFx4vmJTe18mE5zCY52sEAAAAAHCsfBkVKpUSX369HGyBHirHUvswLToUykn3T4XHC417W9PzavW1vS4d7pKAAAABBBHwxNwhNwD0rKJ0uz3pC2zCs8n1LVHmVoMlqJO3FQEAAAAKEmEJzhv9xpp7ihpwcd2MwgjKERqeJEdmmp1ZZQJAAAAjiM8wRn5udKK7+w24+v/KDxfJsluMd7qaim2opMVAgAAAEcgPKF07dtsb2RrNrTN2HnwpEuq38seZarbQwo6sk09AAAA4A0ITyh57nxpzc92A4jVP0ket30+pqLU+hr7KFvd6SoBAACAkyI8oeSk77RHmOaOkVI3FZ6v1c0eZWrYRwoOdbJCAAAAoMgITyheZp8us4bJjDKt+FZy59nnI8pKra6y1zOVr+d0lQAAAMApIzyheGTukRaOs0NTyprC89XaS+2GSY0vlkIjnawQAAAAOCOEJ5zZKNOWOXbHvKVfSnlZ9vmwGKn5FVLboVKlZk5XCQAAABQLwhNOXXa6tOhTac4oaefiwvMVm0ntrpeaXS6FxzpZIQAAAFDsCE8ouh2LpdnvSYs/k3Iy7HMhEVLTAXYDiKpt2MwWAAAAfovwhJPLPWBPyTNrmbbMLjyfUM8OTC0GSVHxTlYIAAAAlArCE45v92p7Wt6Cj6Wsffa5oBCpUV+p7TCpZmdGmQAAABBQCE8olJcjrfzOHmUy7cYLlKkutb1OanmVFFvRyQoBAAAAxxCeIO3dKM0bI837UNqfbJ9zBUn1etlT8+qeJwUFO10lAAAA4CjCU6By50urJ9ujTKt/Mn3H7fMxFaXW10itr5XKJjldJQAAAOA1CE+BJn2HPcJkRppSNxeer93dHmVq0FsKDnWyQgAAAMArEZ4CgdstbfjDHmVa8Z3kzrPPR5aTWl5ph6aEOk5XCQAAAHg1wpM/y9wjLRhrh6Y9awvPJ3WwO+Y1vlgKjXCyQgAAAMBnEJ78jccjbZ5lByazP1N+tn0+LFZqcYXUZqhUqanTVQIAAAA+h/DkL7LSpMWf2nsz7VxSeL5SM3uUqdnlUniMkxUCAAAAPo3w5Ou2L7JHmRZ/JuVk2OdCIqSml9lrmaq2ZjNbAAAAoBgQnnxRTqY9Jc+Epq1zCs+Xr28HphaD7GYQAAAAAIoN4cmX7FolzR0lLfhYykq1zwWFSo36Su2GSTXOZpQJAAAAKCGEJ2+XlyOt+NYeZdrwZ+H5stXt5g+trpJiKjhZIQAAABAQCE/eau8Gae4Yaf6H0v5d9jlXkFT/AntqXp3zpKAgp6sEAAAAAgbhyZu486XVP0mz35PW/Gz6jtvnYypJba6VWl8jlanmdJUAAABAQCI8eYO07fYIkxlpSttSeL72OfYoU4MLpeBQJysEAAAAAh7hySlut7T+d3st04rvJE++fT4yXmp1pb2eKaGO01UCAAAAOIjwVNoy90jzP7K75u1ZV3i++ln2KFOjflJohJMVAgAAADgOwlNp8HikzTPtUaalX0n52fb5sFh7T6a2Q6WKTZyuEgAAAMBJEJ5KUlaatOgTac4oKXlp4fnKLaS2w6SmA6TwGCcrBAAAAFBEhKeSsH2h3TFv8edS7n77XEik1GyAPTWvSms2swUAAAB8DOGpuORkSku/sKfmbZ1beL58Azswmel5kWWdrBAAAADAGSA8naldK+3AtGCclJ1qnwsKlRpfbIemGp0YZQIAAAD8AOHpdORlS8u/sdcybZxaeL5sDbv5Q8urpJhEJysEAAAAUMwIT6di7wZp7mhp3odS5m77nCtIqn+h1O56qfa5UlCQ01UCAAAAKAGEp7+Tnyet/tGemrfmF9N33D4fW1lqfa3U+hqpTFWnqwQAAABQwghPJ5K2zR5hmjdGSttaeL7OufZaJjPaFMz/PgAAACBQ8Or/cG63tH6K3WZ85STJk2+fj4yXWl0ltblOSqjjdJUAAAAAHEB4MvanSAs+shtA7F1feL56J3uUqXE/KSTcyQoBAAAAOMyx7gZLlixRu3btVK5cOd13333yeA6uJSqCffv2qXLlytqwYcPpF2D+vI3TpQk3Si82lCb/2w5O4XFS+5uk22ZI10+Sml9OcAIAAADgTHjKzs5W37591aZNG82ZM0fLli3T6NGji/x8E7Z27Nhxen94Vro0a4T0Vidp1AXS4k+l/Bypckup32vSvSuk3v+TKjQ6va8PAAAAwC85Mm1v0qRJSk1N1YsvvqioqCg99dRTuv322zV06NC/fe4ff/yhr7/+WgkJCaf1Z4e83UEKzjp4J1Jqdpk9Na9q69P6egAA32BmOOS5PcrLd1uTDwAA8InwtHDhQnXs2NEKTkbz5s2t0aeijFjdfPPNevXVV/XAAw/87WPNUSAtLc26deUdkKdiQ7lbXyd3s4FSRBn7Abm5Z/aXgtfLPfg9LrhF4OEa8H+ZOXlavDVNK3aka2NKpjakZGpbapbSDuRq34Fc5eYXpKYQ3TtzsspFhyohOkwJMWGqlRCtuhWiVTcxRk2qxCk2gmXB/orfBeAaQO5pfu9dnlNZbFRM7r33XmVlZemNN944dC4xMVGrVq2y1kCdyGOPPaYFCxZo4sSJqlmzpqZMmWLdHs/jjz+u4cOHH3N+0rvDlZ3YXHK5iulvAwBwSp5bWpvm0tJ9Lut2237JrTP//e6SR1WjpTqxHjUo61GDMh6FsAc6APiNzMxMDRkyxJoNFxcXV+TnOfK2WkhIiMLDj2zCEBERYf0lThSeli9frrffflvz588v0p/x0EMP6Z577jli5CkpKUkt+9502lP+4PvvMEyePFk9e/ZUaGio0+XAAVwD/iE3362pa1I0ccF2TVm9S/uzD24rcVCluHA1q1pGNROirKNauUiVjQpV2chQRYeHKCc3V7/9NkWdOndVeo5HKftzlJyepXW7MrVmV4ZW7kjXln1Z2rJf2rLfpd93SNHhwepeL1G9m1XUOQ0SFRpMkvJl/C4A1wBSUlJO63mOhKf4+Hir297h0tPTFRYWdtzHm8Gxm266SU8++aSqVKlSpD/DhLOjA5phfkD4IQlsXAPgGvBNm1Iy9eGMDfpy/lbtzsg5dL58TLjOa1hBXeqXV5sa5VS5TOTfvmiKCZWqJcSc8DrYkZqlWRv2aPraFP2yfKeS07P13ZId1lE+Jkz9W1fToHZJqp0YU+x/T5QefheAayBwhZ7m992R8GRalI8YMeLQ/fXr11vrk0yoOp5NmzZp6tSpWrx4sdVpr2AkyayVMqNRZsgNAOCfZq3fo/emrtNPy3YeavRg1ild3LKq+rWsouZVyygoqHinYlcqE6F+LapYh9vdVAu27NMPS3ZYwW1Xerbe/WOddfRoVFE3d6uttjXKycV0cADwe46Ep65du1rhZ9SoUVaHPdNtr0ePHgoODrb2cIqNjbU+LlC1alUrYB2uc+fOGj9+vFq2bOnA3wAAUNLmbtyj539cpenrCqdWdK2fqGs61lC3Upw6Z4JZ6+rlrOO+Xg00ZeUujZ+1Sb+uTNbPy3dahxnxuvf8+upUp3yp1AQAcIZja55GjhypwYMHWyNJQUFBVvMHw6x5MuuaDg9F5vFHN4Yw56pVq6aYGKZMAIA/WbwlVS9MXmmFFCMsOEgD2lTT9WfXVL2KsY7WZgJbz8YVrWPtrgyN/HOdJszbqrkb92rIiJnqXLe8/tWrgVomlXW0TgBAyXCsD2u/fv20du1azZ0712pbXtDEoajN/zZs2FDCFQIAStPujGw998MKfTpni3U/OMilgW2r6Y5z66lq2ZOvY3JCncQYPd2/uf7Zo77e+G2Nxs7apKlrdltH/9ZV9eCFDVUhNsLpMgEAxcjRTSwqVaqkPn36OFkCAMBhZtPaj2Zs1AuTVyk9K886d0nLKrq7R33VLB8tb1chLkLDL26qG7rU1ks/r9IX87Zax09Ld+ruHvV0baeadOcDAD/BDoAAAMcs2Zqq+z5fpOXb7Y3Mm1aN0/B+Ta01RL4mKT5KLw5sqWvOqql/T1yiRVtS9eR3y/XZnC167rLmasFUPgDweYQnAECpy8lz6/VfV+uNKWuV7/aoTGSo1YxhcPvq1nQ9X2bWO31129n6dM5mPfvDCq3cma7+b/2lW7rV1j/Oq6fwkMKGSAAA30J4AgCU+mjTvz5bqBU70q37vZtV0hMXN7X2a/IXpkPfoPbVdX6TSnrs66X6ZuE2vfHbWk1etlPPX95CzasxCgUAvohJ2ACAUuF2e6zudJe+Oc0KTvHRYXpjSGu9eWUbvwpOhzN/x9cGt9LbV7W2NtddtTND/d/8S+/8vtb6/wEA8C2MPAEASlxKRrY12vTbwfbj5zeuqKf6N/Pb0HS0C5pWVvtaCXr0qyX6bvF2PT1phaatTdELl7dQYmxg/D8AAH/AyBMAoERNX5uiC1/50wpOYSFB+s8lTfXO1f472nSyUajXh7TS0/2bKTwkSH+s2mX9f5m2ZrfTpQEAiojwBAAoEWbfPjNN76r3Zio5PVt1K8Ro4u1n6+qONeRy+XZTiNNl/t6mKcY3d3ZW/Yox1t5WV78305rGV9R9DgEAziE8AQCK3YGcfN39yQKrVbfppte/VVV9fcfZalQ5zunSvEL9irGaeHtnXd6mmszSJzON745x85WZY+9zBQDwTqx5AgAUq817MnXTh3OtvZtM2/FH+zSyNooN1NGmE4kMC7b2f2qeVFbDv16q7xZt19rkDGtKY40E798cGAACESNPAIBi89fa3er7+lQrOJnuch/f0EHXnV2L4HQC5v+LmcY47qaOVuMI04Ww72tT9edqu7EGAMC7EJ4AAMXi87lbdM17s7QvM1ctqpWx1vV0rJ3gdFk+oV3NeH17Z2e1ql5WaVl5um7UbI2btcnpsgAARyE8AQDOiGl08OJPK61W5Hlujy5qXlmf3HyWKpeJdLo0n1IxLkLjb+qoS1tVtdaJPfTFYj39/XL2gwIAL0J4AgCctuw8uzHEq7+use7ffk4dvTqolSJCg50uzSeFhwTrxYEtdHePetb9d/5Yp9s+nmc14AAAOI/wBAA4LXv35+iqkTM1ccE2hQS59NyA5rqvV0MFBbG+6UzXQd3do75euqKFwoKD9MPSHRr07nRro2EAgLMITwCAU7Zlb6YGvPWXZm/Yq9jwEI0e2l4D2yU5XZZfubRVNX10QweVjQrVwi2puuzt6VYnQwCAcwhPAIBTsmpnuhWc1u3er6plIzXhtk7qXK+802X5pfa14jXh1k7W/+f1u/db/99NJ0MAgDMITwCAIpu7ca8uf3u6dqZlq16FGOuFvdnwFSWnTqL9/7lBxVglp2dr4DvTNXNditNlAUBAIjwBAIpkyspka41T6oFcq6X2Z7ecpUplIpwuKyCY/8+f3nKW2teMV3pWnq5+f5Z+XLrD6bIAIOAQngAAf2vigq26YcwcHcjNV7f6idbmt2WjwpwuK6CUiQzVB8Paq2fjisrJc+vWj+bq09mbnS4LAAIK4QkAcFKjp63XXeMXWHs49WtRRSOuaauosBCnywpIpgX8W1e21qB2STLbP90/YZFGTVvvdFkAEDAITwCAE3rjtzV6/Jtl1sfXdaqpl69oqbAQ/ulwUkhwkJ7u30w3dqll3R/+zTLr+wQAKHm8dQgAOIbH49Erv6zWyz+vtu7fdV49a+NWswcRnGe+D//Xu5Giw0Os79H/flyp/dl5uq9XA75HAFCCePsQAHBMcDIvxguC0/0XNNA/e9bnRbmXbqb7cO9G1v03p6y1RqHcZj4fAKBEEJ4AAEcEp/9+t9x6IW480qeRbute1+mycBI3dq2tJy9pKpNtR/+1QQ9MWKR8AhQAlAjCEwDAYkYsHv96qUZOtRsQPHFxE93QpbbTZaEIrupYQy9c3kJBLumzuVt01/j5ys13O10WAPgdwhMAwApOD3+1WGOmb7RGMJ7p30zXnFXT6bJwCvq3rqY3hrRWaLBL3y7abrUyz87Ld7osAPArhCcACHBmitd9ny/SuFmbrZGL5y9roUHtqztdFk7Dhc0q691r2io8JEg/L0/WjR/M1YEcAhQAFBfCEwAEsLx8t/75yQJNmLdFwUEuvTyolQa0qeZ0WTgD5zSooFHXtVNkaLD+WLVLQ0fPsjrxAQDOHOEJAAJUTp5bd46br68XblNIkEuvD25lbYIL39epbnl9OKy9YsJDNGPdHl3z/iylZeU6XRYA+DzCEwAEILMW5raP52rSkh0KCw7S21e1saZ8wX+0rRmvj27ooLiIEM3duFdXjZypfZk5TpcFAD6N8AQAASYrN183fTDXWhNj1sa8e00b9Whc0emyUAJaJpXVuJs6Kj46TIu2pGrQuzO0OyPb6bIAwGcRngAggGTm5GnYmNn6fdUua02MWRvTvUEFp8tCCWpSpYzG39RRibHhWrEj3QpQO9OynC4LAHwS4QkAAkRGdp6ue3+2pq1JUXRYsMZc395aGwP/V79irD65qaMql4nQmuQMXfHOdG3dd8DpsgDA5xCeACAAmGYB17w3U7M27FFseIg+GNZB7WvFO10WSlHtxBh9evNZqlYuUhtSMjXw7enalJLpdFkA4FMITwDg50yTANMsYN6mfSoTGaqPb+ygNjXKOV0WHJAUH2UFqFrlo62Rp4HvTNfaXRlOlwUAPoPwBAB+bM/+HA0ZMdNqFmCaBoy9sYOaVyvrdFlwUJWykdYUvnoVYrQjLUtXvDNDK3ekO10WAPgEwhMA+Kld6dka/O4MLduepvIx4Rp3Y0ereQBQIS7CaiLRqHKc1X1v0LvTtWRrqtNlAYDXIzwBgB/akZqlK96drpU701UxLlyf3NxRDSrFOl0WvEhCTLjG39hRLaqV0d7MXA0ZMUPzN+11uiwA8GqEJwDwM1v2ZlprWdbt2q+qZSOtNS51EmOcLgteqExUqLWRbtsa5ZSWlaer35ulWev3OF0WAHgtwhMA+JENu/fbXdT2ZKpGQpQ14lQjIdrpsuDFYiNCrbb1Z9VOsNrZX/v+LE1bs9vpsgDAKzkWnpYsWaJ27dqpXLlyuu++++TxeP72OcOHD1d8fLzCw8N16aWXKj2dBa4AUGBNcro14rQtNUu1E6P1yU2mLXWU02XBB0SHh2jU0HbqVj9RB3LzNXT0bP22ItnpsgDA6zgSnrKzs9W3b1+1adNGc+bM0bJlyzR69OiTPufjjz+2jh9++EFLly7V8uXL9cwzz5RazQDgzZZvT7O6piWnZ6uBtSHqWapUJsLpsuBDIkKD9e41bdSzcUXl5Ll104dz9MOSHU6XBQBexZHwNGnSJKWmpurFF19UnTp19NRTT+m999476XM2b96sMWPGqH379qpbt66uuOIKzZ8/v9RqBgBvtXhLqgaPmKGU/TlqUiVO427qqMTYcKfLgg8KDwnWm1e2Vp/mlZWb79HtY+fp64XbnC4LALxGiBN/6MKFC9WxY0dFRdnTSZo3b26NPp3Mgw8+eMT9lStXql69eicd3TJHgbS0NOs2NzfXOhB4Cr7vfP8Dlz9eA/M37dOwD+cpPSvP6pr2/jWtFRvm8qu/Y3Hzx+uguD3fv4lCXdJXC7fr7vHzdSA7R/1bVZW/4BoA1wByT/N77/IUZbFRMbv33nuVlZWlN95449C5xMRErVq1yloD9XfM45o1a6Z58+apSZMmx33M448/bq2ROtrYsWMPhTYA8GVrUqV3VwQr2+1SnViPbmqUr4hgp6uCv3B7pE/XBWl6sj1JZWDtfJ1dsdRfMgBAicjMzNSQIUOs2XBxcXHePfIUEhJiNX04XEREhPWX+Lvw5Ha7df311+uGG244YXAyHnroId1zzz1HjDwlJSXpnHPOUUJCQjH8LeCL7zBMnjxZPXv2VGhoqNPlwAH+dA38siJZ736ySNlutzrViddbQ1oqKsyRX+k+x5+ug5LWx+PRf75fqQ9nbNKn64JVr2EDXXdWDfk6rgFwDSAlJeW0nufIv7SmY57ptnc40zkvLCzsb5/7n//8R3v27NH//ve/kz7OhLOjA5phfkD4IQlsXAPw9Wvgy/lb9K/PFinf7VGPRhX0+pDW1mJ/BNZ1UFqeuLiposJD9M7v6/Tf71cq1y3d1r2u/AHXALgGAlfoaX7fHWkYYVqUT58+/dD99evXW+uTTKg6mW+++cZqMjFhwgSm3gEISKOmrdc/P1loBSezBuWtq9oQnFCiXC6XHrygoe7uYa8zfu6HlXpp8qoibTECAP7GkfDUtWtXaxrdqFGjrPum216PHj0UHBysffv2KT8//5jnmNbkgwcP1muvvWZNv8vIyLCm+QFAIDAvVM0L1uHf2M11rutUU89f3kKhwex1jtIJUHf3qK8HLmho3X/ll9V65ocVBCgAASfIqTVPI0eO1B133KHy5ctr4sSJevbZZ63PmTVPixcvPuY57777rvbv369rr71WsbGx1tG4cWMHqgeA0uV2e6zQZF6wGv/sUV+P9W2soCCX06UhwNzavY7+fZH9b6+ZxmeuS3N9AkCgcGx1cb9+/bR27VrNnTvXalte0MThRO9ivfTSS9YBAIEkN9+t+z9fpC/nb7XuD+/XRNd2qul0WQhg13eupfDQID385RKN/muD9uzP0f8ub27tEQUA/s7R1kyVKlVSnz59nCwBALzW/uw83TF2nn5buUvBQS69cHkLXeJHe+3Ad13ZoYaiwoJ132eLrE10d2dk6+2r2ygugoX3APwbk+UBwAslp2fpinenW8EpIjRI717dhuAEr3Jpq2oaNbSdosOC9dfaFA18e7p2pmU5XRYAlCjCEwB4mTXJGer/5l9asjVN8dFhGndjR53XqKLTZQHH6FIvUZ/cfJYSY8O1Yke6dd2uSU53uiwAKDGEJwDwIrPW79GAt/7Slr0HVDMhSl/c2kmtqp9883DASU2rlrGu09rlo7V13wENeGu6Zm/Y43RZAFAiCE8A4CW+XbRNV42cqdQDuWpVvawm3NpJNctHO10W8LeS4qP0uRX0y1rX75UjZ1rXMwD4G8ITADjMdBl994+1umPsfOXku3V+44oae0NHJcSEO10aUGRmiqm5bns0qqicPLd1Pb/6y2r2ggLgVwhPAOCg7Lx83ff5Ij31/YpDm9++dVUbRYbR9hm+x1y371zdRjd0rmXdf3HyKt39yQJl5eY7XRoA+H6rcgAIZKa9880fztXcjXtl9rt9pE9jDT27plwuNr+F7zJt9R+5qLHqVIjRo18t0cQF27RpT6bevbqt1VgCAHwZI08A4IBl29J08evTrOAUGxGiUUPbW5uPEpzgLwa3r64PhrVXmchQzd+0T5e8MU3Lt6c5XRYAnBHCEwCUsh+W7LA66pnOZLXKR+ur289Wt/qJTpcFFLtOdcrry9s6Wde5ud5NK3OzqS4A+CrCEwCUkny3Ry/8tFK3fDRXB3Lz1blueX1129mqkxjjdGlAiamdGGMFqC71ylvX/T/GzdeT3y5TXr7b6dIA4JQRngCgFKRkZOva92fptV/XHGoMMXpoO5WJCnW6NKDElY0K0+ih7XVr9zrW/ZFT1+uq92Za6/4AwJcQngCghM3btFcXvTZVU9fsVmRosF4Z1FKP92uikGB+BSOwGkk8cEFDvX1Va0WHBWvGuj266NWpmr9pr9OlAUCR8S83AJQQs7/N6GnrdcU707U9NUu1E6M18Y6zdXHLqk6XBjjmgqaVrZ+DOonR2pGWpYHvTNd7U9ezHxQAn0B4AoASkHog19ok9PFvlik336M+zSrr6zs6q37FWKdLAxxXt0Ks1SjlgiaVrJ+P/3y7TMPGzLGmtwKANyM8AUAxm71hj3q/8qe+W7xdIUEuPXpRY70+pJViwtlaDygQGxGqt65qrf9c3ERhIUH6dUWyLnzlT/21drfTpQHACRGeAKCYmO5hL05eZU3TM22ZayRE6fNbO2kY+zcBx2V+Lq4+q6Ym3n626laIUXJ6tq4cOdPqSkk3PgDeiPAEAMVgY8p+XfHuDL36y2q5PdKA1tX03T+6qGVSWadLA7xeo8px+vqOs3VF2ySZpU+mK+WAt6drTXKG06UBwBEITwBwBtxuuynEBS//qbkb9yo2PMTqpvfCwBZM0wNOQVRYiJ69rLleG9xKsREhWrh5n/q8+qdG/rnO2iMNALwB/7IDwBmMNt33+SLNWr/Hut+xdrz+d1kLJcVHOV0a4LP6tqiitjXL6YEJi/XHql168rvl+nHpDj1/eQvVSIh2ujwAAY6RJwA4ReZd8FEHR5tMcIoKC7YWvY+9oSPBCSgGlctEaszQdnq6fzNrT6jZG/ZaP28j/ljHWigAjmLkCQBOweItqXr4q8VatCXVus9oE1ByzSQGt6+uznXL6/7PF2n6uhT99/vl+nL+Vj3VvxnrCQE4gvAEAEWQlpWrF35cqQ9nbLQaQpg1Gfdf0FBXtq+uoCA66QElxbwx8fENHfT53C16atJyLduepkvfnKarO9bQv3o1UFxEqNMlAggghCcA+JuGEBMXbtVT36/QrnR7A8+LW1bRw30aqUJshNPlAQHBvEExsF2SzmtUwRp9+mLeVn0wfaO+X7xD9/dqoAFtqimYNzEAlALCEwCcZLPbJ79dpoUHp+jVKh+t/1zcVJ3rlXe6NCAgJcSE68WBLXVZ62p65KslWrd7v+6fsEhjpm/QY32bqH2teKdLBODnCE8AcJRNKZl65ofl1rvahlmwfts5da3NbiNCg50uDwh4neqW1w93d9UH0zfolV9Wa+m2NA18Z7r6NKus+y9oQFc+ACWG8AQAByWnZemN39Zo3KzNysl3y8wCuqJdku7p2UCJseFOlwfgMGEhQbqhS21d2qqqXpi8SuNnbdJ3i7dbbc3NFL87z61rde0DgOJEeAIQ8HZnZOvtKWutZhDZeXYbZNPhy6xralQ5zunyAPzNVL6nLm1mNZB4ZtIK/b5ql8bO3GQ1mDDnbutex3oMABQHwhOAgJWcnqX3p26wpv5k5uRb59rUKKd7e9a3pgUB8B3mjY4x17e39l57/seVmrVhj96bul7jZm2yWp7f2KW2KpWhyQuAM0N4AhBw1u/er1HTN2nC3K3W9DyjebUyuqdnfXWrn2jtLwPAN5mmEZ/c3FF/rN6tF35aae3JZkKUeZOkf6tqurlbbSWVZSQKwOkhPAEICB6PR3M37tV7K4O0eMY0eTz2+VbVy+q27nXVo1EFQhPgJ8zPsnkjpGu98tY0vjenrLVGpD6Zs1mfzt2sHg0rqIHLZf1eAIBTQXgC4Nf2Z+fpqwVb9eH0jVqxI93sGGOdP69hBd3SvY7a1ihHaAL8lPnZ7t6ggnXM3bhHb01Zq5+XJ2uyORSsn16frmvPrmk1nYgK4yURgL/HbwoAfse8m2xaF5sF4xPmblF6dp51PjwkSC3L5enfV5ytJtXYDwYIJG1qxGvktfFatTNdo6au04S5m7UqOUMPf7nEajRxScuquqxNNWsKL2+oADgRwhMAv7EzLUtfzd+qL+Zt1cqdZpTJVjMhSld1rKGLm1fSX1Mmq37FWEfrBOAc8/P/RL/Gaq4NSi/fRGNnbdaGlEyr26Y56leM0eVtknRxqyqqEEuDCQBHIjwB8Gl79+fo5+U79fXCbZq2ZrfcnsI9YHo2rqiBbZPUpW55BQW5lJub63S5ALxEVIh0WacauqFLHU1bu1ufzdli7RG1ameG/vv9cj09abk61k5Q72aVdUHTSipPu3MAhCcAvmh76gH9tHSnfliyw2pHnF+QmCS1q1lO/VtXs17wlIkMdbROAN7PvLHSpV6idaQeyNW3i7ZZQWrB5n36a22Kdfx74hJ1qJWg3s0rq2ejirQ8BwIY4QmA18vNd2v+pn36c/Uu/bFqlxZuST3i840rx1nvDF/csopqJEQ7VicA32becLmyQw3r2JSSqe+XbNf3i7db7c6nr0uxjke/WqKGlWKtJhTnNEhU6xrlFBpsN6IB4P8ITwC8jtvt0ZpdGZqxLkV/rNpt3WYcbPpgmLXcbaqXswJTryaVlBQf5Wi9APxP9YQo3dKtjnVs3pNphahJS3Zo4ZZ9VudOc7z9+1rFhoeoU90Ea2TK7DFlNusNDqLhBOCvCE8AHJeVm6/FW1M1e8Mezdmw19qPyUyfOVx8dJg61y2vrmbvlvrlWcgNoNSYN2hu7lbHOvbsz7FGwaes3GXtIWXu/7h0p3UYsREhalczXh1qxatlUlk1qVpGMeG83AL8BT/NAEp936Xl29O0ZGuqlmyzb1cnZxyxbsmIDA22NrDtXK+8utZLtKbmmbUJAOAk80bOxS2rWocZJV+0NVV/rd1tbcJr3vxJz8rTryuSraNgpLxuYoyaVSujFtXKWrdm2h/7SgG+iZ9cACUiLStX63bt19rkDGsKXsHt+t375TkyJ1lMJyvT7KFtzXjr1kx9YR0BAG9m3tAxo0vmuK27lJfv1rLtaVaQMocZUd+emmW9QWQOs41CQaCqVi5S9SvEql7FWKs9ummhXjsxmlAFeDl+QgGcFvMiwbwo2LL3gLbuO6AtezOtj82tCU3J6dknfG7FuHA1rVLGms7StEqcmlYto8plItiYEoBPCwkOUvNqZa3jhi61rXPJ6VlavCXVanSzeMs+K1DtzsjR5j0HrOOXgyNUBcrHhKl6fJR9JEQf+rhquUglxoRb2zAACMDwtGTJEg0dOlRr1qzRDTfcoOeee+5vXzh9/vnnuvfee629Wl544QUNHjy41OoFAoGZOpeelat9mblK2Z+jXenZ2pWRbd9aR5Z1m3zwOHqq3dEqxIarTmKM6laIUZ3EaNWpEKMGlWJZrwQgYJjfd+c1MkfFQ+dSMrKt/aRWJ6dr1U5zZGj1znTtzcy1gpU55m3ad9yvZ8KV+ZqmXbp5I6rgYzOdsFxUmOKjQ1U2KkxlI0OtMAfAD8JTdna2+vbtq169emn8+PH6xz/+odGjR1th6mRh68orr9Qbb7yhDh06qH///mrdurUaNGhQqrUD3jgClJPvVlau21pPdCA337rNzMk/eNgfW58ztzn51pS61MxcqynDvgM51q25n5ZV2NGuKMw7oFXLRlrTT8xhPjbvjtZMsINSXAT7LAHA0RJiwnWWOeokHHHe/B7evDdTG1MytWmPOfZbt+b+zrQs5eZ7DoUrMz3w78RFhFihygpTUaFW4wpzRB+8PfRxhPk4WNFh9v2I0CCFhwQrvOA2xNwGMTsAcCo8TZo0SampqXrxxRcVFRWlp556SrfffvtJw9PIkSN1zjnnWKNUxh133KEPP/xQTz755Cn92T8vT1ZcmZxD94/3vvnx1mN4jvvIEz32eI87+Tv0p/Pnn+hLnmlNxz17pjUd988v2tc88WOL+MCDj83Pz9fSHS7tmblJwcHBx/+axzl5ou+cGXQxi4XzPR5rBKbg48Jz9tfLP875Yx/rkfvgbV6+xwpDOXnuwts8t7XXUcG57MPu/83gz2mJCgu23sGsEBduTRNJjLXf3TS3BUeVMhHWOiWaOABA8SgTFaoyUWWsqcxHM/9e7M3M0Y60LCWnZVu3Ow8d2VbXv32ZOdZtwRth5tb6OCWzWOorCFERoYXBygStkCBzuBQS7LI+Nq3aQ4Nd1q31uYMfh5rPWY858rz5VyTI5ZL558QENJPRjne/4HGug+fN54++bx5VcN484eh/oQoCYH5+nhbtcil7/jbrNYH9ucMfd9TzjvpKJ8uRh4fMY//8E3/dYz9XtD8Ppydt3/FHd70yPC1cuFAdO3a0gpPRvHlzLVu27G+fc+GFFx663759ez3xxBMnHd0yR4G0NPsdmvsmLFFQOHvCBK5gfb5+hfyV+UfMdKmLDgtW5MEjOizEOmcCUdTBc2UiQhUXGWJN64iLDD10WyYyxBotKuqcevOPT36+fIaZ8nv4LQIT1wF89RqICw9SXGKU6idG/e2MBDOjwEwD3Gdu9+cqNSvXmoGQkW3PRNifk6eMLDMbIe+I82amgnmDLjsvX1l57iPeaLTPu095loL3CtbHa5Y4XQQc4s7O9J3wZIJMrVq1jkjoJvXv3btX5cqVK9Jz4uLitG3bthP+GU8//bSGDx9+zPmaMR6FRBz5Vn1Rw/yJU/+xb/2fyhsERX034XgPK4k/51S+7vG+5inVVEr/n054/gzqN4+z3w2TTNSw3yE77OMTft5z7GMP+zjYJZnsEmI+Pnhr3/dYt4d//vDb0CD7+UWSdfDYKx2QfWxX4Jg8ebLTJcALcB0gkK6BqINHYsEJM9gSffA4AROc8j1SnlvKcUt5HinX3Lrt21y3y7p1H3zc0bcFH7utNbWH3XpcMu+7ud2ybs3LKPNKyjqO9/Fh58x/zNcpCHUnfM5xJqOcbKLGySYIHf2p4vo6h7/iOLXnoTjkhnq02VfCU0hIiMLDw484FxERoczMzBOGp6OfU/D4E3nooYd0zz33HBG+kpKS9OntXZWQcOQcYwQG8w6j+YeyZ8+eCg1lLU4g4hqAwXUArgFwDSAlJUWVT231j3PhKT4+3moAcbj09HSFhYWd9Dm7du0q8uNN0Do6oBnmB4QfksDGNQCuARhcB+AaANdA4Ao9ze+7Iz0s27Vrp+nTpx+6v379emt9kglIRX3O/PnzVbVq1RKvFQAAAAAcC09du3a1ptGNGjXKum+67fXo0cNa97Rv3z6rK9rRBgwYYLU1X7x4sTIyMvTqq69arc4BAAAAwG/Dk1m/ZFqPm3bj5cuX18SJE/Xss89anzNrnkxAOlqLFi101113qW3bttaIkwlat912mwPVAwAAAAhEjqx5Mvr166e1a9dq7ty5VtvygiYOJ9sP6b///a+1Ue7WrVvVrVu3k655AgAAAAC/CE9GpUqV1KdPn1N6TuPGja0DAAAAAPx+2h4AAAAA+BrCEwAAAAAUAeEJAAAAAIqA8AQAAAAARUB4AgAAAIAiIDwBAAAAQBEQngAAAACgCAhPAAAAAFAEhCcAAAAAKALCEwAAAAAUAeEJAAAAAIqA8AQAAAAARRCiAOHxeKzb9PR0hYaGOl0OHJCbm6vMzEylpaVxDQQorgEYXAfgGgDXANLT04/ICEUVMOEpJSXFuq1Vq5bTpQAAAADwkoxQpkyZIj8+YMJTfHy8dbtp06ZT+h8E/2HeXUpKStLmzZsVFxfndDlwANcADK4DcA2AawCpqamqXr36oYxQVAETnoKC7OVdJjjxQxLYzPefayCwcQ3A4DoA1wC4BhB0MCMU+fElVgkAAAAA+BHCEwAAAAAUQcCEp/DwcD322GPWLQIT1wC4BmBwHYBrAFwDCD/Na8DlOdX+fAAAAAAQgAJm5AkAAAAAzgThCQAAAACKgPAEAAAAAEUQkOHpgQceUN++fZ0uA17gggsu0OjRo50uA6Vs4sSJql27tkJCQtSyZUstX77c6ZIAlBJ+/nE4XgcEtgdOIxMEXHhatGiR3nzzTb3yyitOlwKHffzxx/rxxx+dLgOlbO3atRo6dKieeeYZbd26VfXr19cNN9zgdFkoJUuWLFG7du1Urlw53XfffaJnUmDh5x+H43VAYFt0mpkgoMKT2+3WTTfdpH/+85/Wu04IXHv27NG9996rBg0aOF0KSpl5l9m8cBo4cKAqVqyoW2+9VfPnz3e6LJSC7Oxs6x3GNm3aaM6cOVq2bBnvOAcYfv5RgNcBgc19BpkgoMLT22+/rcWLF6tmzZr6+uuvlZOT43RJcIj5hXnppZeqY8eOTpeCUnbRRRdZvzALrFy5UvXq1XO0JpSOSZMmKTU1VS+++KLq1Kmjp556Su+9957TZaEU8fOPArwOCGxvn0Em8LvwdMkll6hs2bLHHK+//rq1EZZJlxs3btRLL72kzp0768CBA06XjFK+Dn777Tf98ssveu6555wuEw5dAwXML8sXXnhBt9xyi6O1onQsXLjQeqEUFRVl3W/evLk1+oTAxM9/4OJ1QGDLyMg4o0wQIj/zzjvvHPcv/8cff2j//v3WD0z58uWVl5enZs2a6cMPPzziXSj493UQHx+vtm3b6q233lJsbKwjtcH5a6CA+eUZHR3NmocAkZaWplq1ah2673K5FBwcrL1791proBBY+PkPTFlZWbr55pt5HRDAvvjiizPKBH4Xnswc5uMZO3as9Y6j+Z9kmC475l3HNWvWlHKFcPI6ePjhh63F4n369Cn1muAd10CBX3/9VW+88YZmzJih0NDQUqsLzjG/98PDw484FxERoczMTMJTgOHnP3D95z//4XVAgNuyZcsZZQK/C08nUq1atWPehTZDdZ06dXKsJpQ+E6J37dplTd8yzIumTz/9VLNmzbI6riAwrF+/XoMHD7ZePDVu3NjpclBKzKij6bZ3uPT0dIWFhTlWE0ofP/+BjdcBqHaGmcDlCZA+rSkpKdbcxmeffdZaMGqG7Exv91WrVikpKcnp8lCK7zaY4dkC//rXv6x3H6677rpD70DAv5lfmGbq5tlnn201Dihgpu+YaVzw79EGMyWj4N1F8yLavHg289/N9D34P37+wesApJxhJvC7hhEnkpCQoO+//15jxoyx9nUwPd3NOw0Ep8B7t8F0Vik4YmJirF+W/MIMHD/99JPVJGDEiBHWfPeCw7zrBP/WtWtXa93TqFGjrPum216PHj0ITgGEn3/wOgAJZ5gJAmbkCQAA05LWTNmKjIxUUFCQpkyZwtQtAECREZ4AAAFlx44dmjt3rjVVx7wDCQBAURGeAAAAAKAIAmbNEwAAAACcCcITAAAAABQB4QkAAAAAioDwBAAAAABFQHgCAAAAgCIgPAEAAABAERCeAAAAAKAICE8AAAAAUASEJwCA39u+fbvKlCmjGTNmWPcfe+wxnX322WKfeADAqXB5+JcDABAAXnzxRU2cOFGfffaZ6tevrylTpqhly5ZOlwUA8CGEJwBAQMjLy7PCUlRUlDp06KDXXnvN6ZIAAD6GaXsAgIAQEhKiG264QbNnz9Ytt9zidDkAAB/EyBMAICCkpaWpcePGatKkiWJiYjRhwgSnSwIA+BhGngAAAeHhhx9Wp06drNBkGkd8++23TpcEAPAxjDwBAPzerFmzdO6552rp0qWqUaOGPv74YytMmfvR0dFOlwcA8BGEJwAAAAAoAqbtAQAAAEAREJ4AAAAAoAgITwAAAABQBIQnAAAAACgCwhMAAAAAFAHhCQAAAACKgPAEAAAAAEVAeAIAAACAIiA8AQAAAEAREJ4AAAAAQH/v/wFg2mJv8sB0CwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x600 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# PAO分布\n",
    "from scipy._lib._util import _lazyselect, _lazywhere\n",
    "\n",
    "_norm_pdf_C = np.sqrt(2*np.pi)\n",
    "_norm_pdf_logC = np.log(_norm_pdf_C)\n",
    "\n",
    "def _norm_pdf(x):\n",
    "\treturn np.exp(-x**2/2.0) / _norm_pdf_C\n",
    "\n",
    "def _norm_cdf(x):\n",
    "\treturn sc.ndtr(x)\t\n",
    "\n",
    "def pdf_pao1(x, sigma: float, skew: float, hurst: float) -> None:\n",
    "\t# skew: 偏度, 0≤skew<2\n",
    "\t# hurst: 峰度, ∞>hurst≥1.772\n",
    "\tc\n",
    "\t# (a+1)/4 : 确保a==3时，分布为正态分布; a==0时, 分布为均匀分布且y不等于0\n",
    "\treturn (a+1)/4 * (np.sqrt(2*np.pi)*sigma) * np.exp(-a*(x/(2*sigma))**2)\n",
    "\n",
    "def pdf_pao(x, sigma: float, skew: float, hurst: float) -> None:\n",
    "\t# skew: 偏度, 0≤skew<2\n",
    "\t# hurst: 峰度, ∞>hurst≥1.772\n",
    "\t# (a+1)/4 : 确保a==3时，分布为正态分布; a==0时, 分布为均匀分布且y不等于0\n",
    "\treturn 1/np.sqrt(2*np.pi) * np.exp(-x**2/2)\n",
    "\n",
    "x = np.linspace(-6, 6, 1000)\n",
    "\n",
    "def plt_pao(x, sigma: float, skew: float, hurst: float) -> None:\n",
    "\tplt.plot(x, pdf_pao(x, sigma, skew, hurst), label=f'sigma={sigma:.3f}, skew={skew:.3f}, hurst={hurst:.3f}')\n",
    "\n",
    "plt.figure(figsize=(10, 6))\n",
    "# plt_pao(x, 1, 1, 1.772)\n",
    "# plt_pao(x, 1, 3, 3)\n",
    "# plt_pao(x, 1, 1, 4)\n",
    "skew = 0.005\n",
    "hurst = 3\n",
    "a = 1 / 2\n",
    "b = skew\n",
    "c = (hurst-3)\n",
    "d = 1\n",
    "\n",
    "b = -0.1\n",
    "c = 1\n",
    "plt.plot(x, 1/np.sqrt(2*np.pi) * np.exp(-x**2/2))\n",
    "plt.plot(x, 1/np.sqrt(2*np.pi) * np.exp(-a*x**2) + x/10 + 1)\n",
    "\n",
    "plt.legend()\n",
    "plt.grid()\n",
    "plt.title('PAO分布')\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('f(x)')\n",
    "plt.xlim(-6, 6)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.13.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
